Activity Budget Analysis
Cali Wilson
5/20/2022
How does urbanization influence ibis behaviors relevant to transmission?
Specifically, how does variation in flock size and frequency of provisioning influence these behaviors?
Data:
10-minute focal follow observations of birds at 5 urban parks in South Florida. Follows conducted during summer 2019.
Focal follows were only conducted in urban parks
Note: Focal follows were only conducted when people were not actively feeding the birds so foraging means probing on the ground for food. The few instances of consuming human provided food are cases where ibis were foraging ‘naturally’ and happened upon scraps of anthropogenic food that was already on the ground.
Sample Sizes / Data Structure
Observations consist of 168 focal follows lasting > 2 minutes. Five urban parks were samples for ~30 focals per site. Sites were visited 2-4 times each.
Provisioning estimates: Each visit, # of people and groups feeding birds were recorded for duration of visit. Provisioning was estimated by counting the number of groups that fed the birds per visit and dividing by total time I was at the site to get # provisioning events/ hour. Each focal follow was paired with the average # of provisioning events per hour on day/site of the focal follow.
Flock size estimates: Every ~30 minutes, flock size was recorded. I paired each focal follow with the nearest estimate of flock size.
Birds seem to spend the most time foraging, grooming, and being vigilant.
I am going to re-group these variables according to their relevance to different modes of parasite transmission (see different sections below).
Summary Statistics
How many focal follows lasted > 2 minutes? *important- when we add in predictors later, we lose data for 3 follows because they lack corresponding predictor data.
## [1] 150What was the total amount of time spent observing ibis? (this only includes obs >2 minutes aka the ones used in the analysis)
Total obs time in seconds:
## [1] 79994Total obs time in minutes:
## [1] 1333.233Mean obs time:
## [1] 533.2933How many different behaviors were observed on average for each ibis during the observation?
## focalID n_behaviors
## Length:150 Min. :2.000
## Class :character 1st Qu.:3.000
## Mode :character Median :4.000
## Mean :3.853
## 3rd Qu.:5.000
## Max. :6.000How many birds showed 2 or more behaviors?
nrow(sum3.1[sum3.1$n_behaviors >= 2,]) #all birds observed for at least 2 minutes did 2 or more behaviors ## [1] 150nrow(sum3.1[sum3.1$n_behaviors < 2,])## [1] 0Most commonly observed behaviors observed?
How many birds foraged?
sum4 <- simp2
sum4f <- subset(sum4, Activity=="Foraging" & activity_count > 0)
summary(sum4f)## focalID focalduration Activity activity_count
## Length:110 Min. :139.0 Length:110 Min. : 4.0
## Class :character 1st Qu.:592.0 Class :character 1st Qu.: 53.5
## Mode :character Median :596.0 Mode :character Median :148.0
## Mean :537.1 Mean :212.6
## 3rd Qu.:597.0 3rd Qu.:362.0
## Max. :598.0 Max. :585.0How many birds groomed?
sum4g <- subset(sum4, Activity=="Grooming" & activity_count > 0)
summary(sum4g)## focalID focalduration Activity activity_count
## Length:113 Min. :122.0 Length:113 Min. : 2.0
## Class :character 1st Qu.:557.0 Class :character 1st Qu.: 22.0
## Mode :character Median :596.0 Mode :character Median :140.0
## Mean :536.5 Mean :170.8
## 3rd Qu.:597.0 3rd Qu.:275.0
## Max. :598.0 Max. :572.0How many birds were vigilant?
sum4v <- subset(sum4, Activity=="Vigilant" & activity_count > 0)
summary(sum4v)## focalID focalduration Activity activity_count
## Length:140 Min. :122.0 Length:140 Min. : 4.00
## Class :character 1st Qu.:579.8 Class :character 1st Qu.: 41.75
## Mode :character Median :596.0 Mode :character Median :132.50
## Mean :539.2 Mean :178.65
## 3rd Qu.:597.0 3rd Qu.:269.00
## Max. :598.0 Max. :591.00What behaviors did ibis spend the most time doing?
sum5<-simp2
sum5 <- sum5 %>% group_by(Activity) %>%dplyr::summarise(sum(activity_count))
sum5## # A tibble: 7 × 2
## Activity `sum(activity_count)`
## <chr> <dbl>
## 1 Drinking 243
## 2 Foraging 23387
## 3 Grooming 19305
## 4 Other 812
## 5 Resting 2934
## 6 Vigilant 25011
## 7 Walking 8302Which behaviors were observed most frequently in ibis?
sum6<-simp2
sum6<- as.data.frame(subset(simp2, activity_count > 0))
sum6.1 <- sum6 %>% group_by(Activity) %>%
dplyr::summarise(n_focals=n(),
.groups = 'drop')
sum6.1## # A tibble: 7 × 2
## Activity n_focals
## <chr> <int>
## 1 Drinking 19
## 2 Foraging 110
## 3 Grooming 113
## 4 Other 61
## 5 Resting 20
## 6 Vigilant 140
## 7 Walking 115Visualizations of predictors
Predictors:
- nearest flock size estimate
- feeding events per hour
- site? (not really interested in site-level differences per say but might help us understand if there are other inherent differences across sites that we aren’t currently measuring)
What is the distribution of provisioning frequency? What about average flock size? Do we see any relationship between # feeding events and average flock size?
Note: Here I’m showing avg flock size and nearest flock size
# Provisioning and Flock Density
ggplot(predictors, aes(x=feedingperhour, y=flocksize_nearest)) +
geom_point() ## Warning: Removed 3 rows containing missing values (`geom_point()`).
ggplot(predictors, aes(x=feedingperhour, y=log(flocksize_nearest))) +
geom_point() ## Warning: Removed 3 rows containing missing values (`geom_point()`).
lm1 <- lm(predictors$flocksize_nearest~predictors$feedingperhour)
summary(lm1)##
## Call:
## lm(formula = predictors$flocksize_nearest ~ predictors$feedingperhour)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.368 -31.945 -16.797 5.128 169.065
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.111 7.049 5.69 6.58e-08 ***
## predictors$feedingperhour 9.647 8.178 1.18 0.24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 53.47 on 148 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.009315, Adjusted R-squared: 0.002621
## F-statistic: 1.392 on 1 and 148 DF, p-value: 0.24Model visualizations and Model Selection
Standardize predictor variables
#library(mosaic) # standardizing variables
ForagePreen_wide <-
ForagePreen_wide %>%
mutate(avg_flock_size_s = scale(avg_flock_size))%>%
mutate(flocksize_nearest_s = scale(flocksize_nearest))%>%
mutate(feedingperhour_s = scale(feedingperhour))
Vig_wide <-
Vig_wide %>%
mutate(avg_flock_size_s = scale(avg_flock_size))%>%
mutate(flocksize_nearest_s = scale(flocksize_nearest))%>%
mutate(feedingperhour_s = scale(feedingperhour))Foraging models:
Main Q:
Does flock size, provisioning, both, or neither influence foraging time?
Foraging includes foraging on soil, in water, on human objects, and consuming human provided food.
Note: Focal follows were only conducted when people were not actively feeding the birds so foraging means probing on the ground for food. The few instances of consuming human provided food are cases where ibis were foraging ‘naturally’ and happened upon scraps of anthropogenic food that was already on the ground.
Exploratory Visualizations
Model selection:
Count data: Does poisson or negative binomial fit better? Build global model (two predictors + interactions) and compare AIC.
pm1 = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s + offset(log(duration)), ForagePreen_wide, family=poisson)
nbm1 = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), ForagePreen_wide, family=nbinom2)
zipm1 = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~flocksize_nearest_s*feedingperhour_s, ForagePreen_wide, family=poisson)
zipm1.int = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~1, ForagePreen_wide, family=poisson)
zinbm1 = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~flocksize_nearest_s*feedingperhour_s, ForagePreen_wide, family=nbinom2)
zinbm1.int = glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~1, ForagePreen_wide, family=nbinom2) #intercept only zero inflated term (equally applies zero inflated to whole model)
AICtab(pm1,nbm1,zipm1,zipm1.int,zinbm1,zinbm1.int) #zinbm1.int is best fit ## dAIC df
## zinbm1.int 0.0 6
## zinbm1 2.5 9
## nbm1 59.9 5
## zipm1.int 14290.9 5
## zipm1 14293.2 8
## pm1 27285.8 4Is this the case for other candidate models (not just the global model)?
#just flock size
pm2 = glmmTMB(Foraging~flocksize_nearest_s + offset(log(duration)), ForagePreen_wide, family=poisson)
nbm2 = glmmTMB(Foraging~flocksize_nearest_s + offset(log(duration)), ForagePreen_wide, family=nbinom2)
zipm2 = glmmTMB(Foraging~flocksize_nearest_s + offset(log(duration)), zi=~flocksize_nearest_s, ForagePreen_wide, family=poisson)
zinbm2 = glmmTMB(Foraging~flocksize_nearest_s + offset(log(duration)), zi=~flocksize_nearest_s, ForagePreen_wide, family=nbinom2)
zinbm2.int = glmmTMB(Foraging~flocksize_nearest_s + offset(log(duration)), zi=~1, ForagePreen_wide, family=nbinom2)
AICtab(pm2,nbm2,zipm2,zinbm2,zinbm2.int) #the two zinbm models are the best ## dAIC df
## zinbm2.int 0.0 4
## zinbm2 1.8 5
## nbm2 58.6 3
## zipm2 15518.7 4
## pm2 29415.6 2#just provisioning
pm3 = glmmTMB(Foraging~feedingperhour_s + offset(log(duration)), ForagePreen_wide, family=poisson)
nbm3 = glmmTMB(Foraging~feedingperhour_s + offset(log(duration)), ForagePreen_wide, family=nbinom2)
zipm3 = glmmTMB(Foraging~feedingperhour_s + offset(log(duration)), zi=~feedingperhour_s, ForagePreen_wide, family=poisson)
zinbm3 = glmmTMB(Foraging~feedingperhour_s + offset(log(duration)), zi=~feedingperhour_s, ForagePreen_wide, family=nbinom2)
zinbm3.int = glmmTMB(Foraging~feedingperhour_s + offset(log(duration)), zi=~1, ForagePreen_wide, family=nbinom2)
AICtab(pm3,nbm3,zipm3,zinbm3,zinbm3.int) #the two zinbm models are the best## dAIC df
## zinbm3.int 0.0 4
## zinbm3 0.2 5
## nbm3 57.8 3
## zipm3 14817.4 4
## pm3 28058.0 2#Flock size and provisioning
pm4 = glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s + offset(log(duration)), ForagePreen_wide, family=poisson)
nbm4 = glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s +offset(log(duration)), ForagePreen_wide, family=nbinom2)
zipm4 = glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s +offset(log(duration)), zi=~flocksize_nearest_s+feedingperhour_s, ForagePreen_wide, family=poisson)
zinbm4 = glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s +offset(log(duration)), zi=~flocksize_nearest_s+feedingperhour_s, ForagePreen_wide, family=nbinom2)
zinbm4.int = glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s +offset(log(duration)), zi=~1, ForagePreen_wide, family=nbinom2) #intercept only zero inflated term
AICtab(pm4,nbm4,zipm4,zinbm4,zinbm4.int) #the two zinbm models are the best ## dAIC df
## zinbm4.int 0.0 5
## zinbm4 2.0 7
## nbm4 60.2 4
## zipm4 14327.2 6
## pm4 27600.9 3Yes- the zero inflated models always outperform the others.
I can also double check dispersion and zero inflation in my models using Dharma:
#check fit of poisson model:
simulationOutput.pm1 <- simulateResiduals(fittedModel = pm1)
plot(simulationOutput.pm1)
#overdispersed or underdispersed?
testDispersion(simulationOutput.pm1,alternative = "greater") #test for overdispersion 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 186.35, p-value < 2.2e-16
## alternative hypothesis: greatertestDispersion(simulationOutput.pm1,alternative = "less") #test for underdispersion 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 186.35, p-value = 1
## alternative hypothesis: less#poisson model is overdispersed
#zero inflated?
testZeroInflation(simulationOutput.pm1) #yes zero inflated 
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = Inf, p-value < 2.2e-16
## alternative hypothesis: two.sided#Poisson model is overdispersed and zero inflated- will the neg binom model do better?
simulationOutput.nbm <- simulateResiduals(fittedModel = nbm1)
plot(simulationOutput.nbm)
#overdispersed or underdispersed?
testDispersion(simulationOutput.nbm,alternative = "greater") #test for overdispersion 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.29045, p-value = 1
## alternative hypothesis: greatertestDispersion(simulationOutput.nbm,alternative = "less") #test for underdispersion 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.29045, p-value < 2.2e-16
## alternative hypothesis: less#This model is underdispersed
testZeroInflation(simulationOutput.nbm) #still slightly zero inflated 
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.5642, p-value = 0.016
## alternative hypothesis: two.sided#Now test a zero inflated neg binom model (look at zinbm1 and zinbm1.int)
simulationOutput.zinbm1 <- simulateResiduals(fittedModel = zinbm1)
plot(simulationOutput.zinbm1)
simulationOutput.zinbm1.int <- simulateResiduals(fittedModel = zinbm1.int)
plot(simulationOutput.zinbm1.int)
testDispersion(simulationOutput.zinbm1,alternative = "greater") #test for overdispersion 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.6559, p-value = 0.924
## alternative hypothesis: greatertestDispersion(simulationOutput.zinbm1,alternative = "less") #test for underdispersion 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.6559, p-value = 0.076
## alternative hypothesis: lesstestDispersion(simulationOutput.zinbm1.int,alternative = "greater") #test for overdispersion 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.67235, p-value = 0.9
## alternative hypothesis: greatertestDispersion(simulationOutput.zinbm1.int,alternative = "less") #test for underdispersion 
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.67235, p-value = 0.1
## alternative hypothesis: less#disperson is fine for both of these, what about zero-inflation?
testZeroInflation(simulationOutput.zinbm1)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0068, p-value = 1
## alternative hypothesis: two.sidedtestZeroInflation(simulationOutput.zinbm1.int)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0021, p-value = 1
## alternative hypothesis: two.sided#Both of these models are good- no dispersion or zero inflation issuesThese results suggest that the zero inflated models make the most sense (they fix the dispersion and zero inflation issues). ZI models with zi=1 and zi= conditional fixed effects are both sufficient in fixing these issues.
For these data, I would think most zeros are false zeros that are due to experimental design AKA 0 foraging does not mean the bird never forages but rather it was not observed for long enough to observe foraging. I don’t think these false zeros could be explained by my predictors. Because of this, I think I should apply a single zero inflation parameter to my entire model (ziformula~1).
Now I will move forward with zero inflated negative binomial models with a single zero inflated term applied to all models.
First build models without any random effects:
#intercept only
fmod0.3 <-glmmTMB(Foraging~1+ offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod0.3)## Family: nbinom2 ( log )
## Formula: Foraging ~ 1 + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1571.0 1580.1 -782.5 1565.0 147
##
##
## Dispersion parameter for nbinom2 family (): 0.987
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.91841 0.09627 -9.54 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0340 0.1881 -5.497 3.85e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning
fmod11 <-glmmTMB(Foraging~feedingperhour_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod11)## Family: nbinom2 ( log )
## Formula: Foraging ~ feedingperhour_s + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1567.9 1579.9 -779.9 1559.9 146
##
##
## Dispersion parameter for nbinom2 family (): 1.03
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.95398 0.09454 -10.091 <2e-16 ***
## feedingperhour_s -0.19827 0.08581 -2.311 0.0209 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0324 0.1878 -5.497 3.86e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#flock size
fmod12 <-glmmTMB(Foraging~flocksize_nearest_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod12)## Family: nbinom2 ( log )
## Formula: Foraging ~ flocksize_nearest_s + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1567.6 1579.6 -779.8 1559.6 146
##
##
## Dispersion parameter for nbinom2 family (): 1.03
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.94968 0.09438 -10.06 <2e-16 ***
## flocksize_nearest_s -0.26320 0.10240 -2.57 0.0102 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0343 0.1881 -5.5 3.8e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning + flock size
fmod13 <-glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod13)## Family: nbinom2 ( log )
## Formula:
## Foraging ~ flocksize_nearest_s + feedingperhour_s + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1565.2 1580.3 -777.6 1555.2 145
##
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.97975 0.09296 -10.539 <2e-16 ***
## flocksize_nearest_s -0.24333 0.10301 -2.362 0.0182 *
## feedingperhour_s -0.18472 0.08670 -2.131 0.0331 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0330 0.1878 -5.499 3.82e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fmod13.1 <-glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod13.1)## Family: nbinom2 ( log )
## Formula:
## Foraging ~ flocksize_nearest_s + feedingperhour_s + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1565.2 1580.3 -777.6 1555.2 145
##
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.97975 0.09296 -10.539 <2e-16 ***
## flocksize_nearest_s -0.24333 0.10301 -2.362 0.0182 *
## feedingperhour_s -0.18472 0.08670 -2.131 0.0331 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0330 0.1878 -5.499 3.82e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning x flock size
fmod14 <-glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod14)## Family: nbinom2 ( log )
## Formula:
## Foraging ~ flocksize_nearest_s * feedingperhour_s + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1567.1 1585.2 -777.6 1555.1 144
##
##
## Dispersion parameter for nbinom2 family (): 1.07
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.97711 0.09345 -10.456 <2e-16 ***
## flocksize_nearest_s -0.24949 0.10501 -2.376 0.0175 *
## feedingperhour_s -0.20032 0.09952 -2.013 0.0441 *
## flocksize_nearest_s:feedingperhour_s -0.05511 0.17401 -0.317 0.7515
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0333 0.1879 -5.499 3.81e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AICtab(fmod0.3,fmod11,fmod12,fmod13,fmod14)## dAIC df
## fmod13 0.0 5
## fmod14 1.9 6
## fmod12 2.4 4
## fmod11 2.7 4
## fmod0.3 5.8 3fmod13 (flock size and feeding) and fmod14 (flock size * feeding) are best fitting models.
Now build the same models but include site as a random effect and see if same models are best fitting. I think site makes sense biologically but it only has 4 levels so questionable as a RE.
#intercept only
fmod0.2 <-glmmTMB(Foraging~1+ (1|site) + offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod0.2)## Family: nbinom2 ( log )
## Formula: Foraging ~ 1 + (1 | site) + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1570.3 1582.4 -781.2 1562.3 146
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.08102 0.2846
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.05
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.9674 0.1599 -6.052 1.43e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.033 0.188 -5.498 3.84e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning
fmod15 <-glmmTMB(Foraging~feedingperhour_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod15)## Family: nbinom2 ( log )
## Formula:
## Foraging ~ feedingperhour_s + (1 | site) + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1569.8 1584.9 -779.9 1559.8 145
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01202 0.1097
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.04
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.9640 0.1126 -8.561 <2e-16 ***
## feedingperhour_s -0.1861 0.1062 -1.752 0.0798 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0327 0.1879 -5.497 3.86e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#flock size
fmod16 <-glmmTMB(Foraging~flocksize_nearest_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod16)## Family: nbinom2 ( log )
## Formula:
## Foraging ~ flocksize_nearest_s + (1 | site) + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1569.6 1584.6 -779.8 1559.6 145
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 4.094e-07 0.0006398
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.03
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.94967 0.09439 -10.062 <2e-16 ***
## flocksize_nearest_s -0.26319 0.10244 -2.569 0.0102 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0343 0.1881 -5.5 3.8e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning + flock size
fmod17 <-glmmTMB(Foraging~flocksize_nearest_s+feedingperhour_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod17)## Family: nbinom2 ( log )
## Formula:
## Foraging ~ flocksize_nearest_s + feedingperhour_s + (1 | site) +
## offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1567.2 1585.3 -777.6 1555.2 144
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 2.283e-09 4.779e-05
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.97975 0.09296 -10.539 <2e-16 ***
## flocksize_nearest_s -0.24333 0.10301 -2.362 0.0182 *
## feedingperhour_s -0.18472 0.08670 -2.131 0.0331 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0330 0.1878 -5.499 3.82e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning x flock size
fmod18 <-glmmTMB(Foraging~flocksize_nearest_s*feedingperhour_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod18)## Family: nbinom2 ( log )
## Formula:
## Foraging ~ flocksize_nearest_s * feedingperhour_s + (1 | site) +
## offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1569.1 1590.2 -777.6 1555.1 143
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 2.117e-09 4.601e-05
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.07
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.97710 0.09345 -10.455 <2e-16 ***
## flocksize_nearest_s -0.24949 0.10501 -2.376 0.0175 *
## feedingperhour_s -0.20032 0.09952 -2.013 0.0441 *
## flocksize_nearest_s:feedingperhour_s -0.05511 0.17401 -0.317 0.7515
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0333 0.1879 -5.499 3.81e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AICtab(fmod0.2,fmod15,fmod16,fmod17,fmod18) #fmod17 and fmod18 = best fit## dAIC df
## fmod17 0.0 6
## fmod18 1.9 7
## fmod16 2.4 5
## fmod15 2.6 5
## fmod0.2 3.1 4Models with flock size and feeding or flock size * feeding are still the best fitting models (fmod17 and fmod18 here).
Summary estimates aren’t too different for the top performing models regardless of whether they include site as a random effect or not.
Model diagnostics
Plotting predicted values:
Plots made using ggeffects: https://strengejacke.github.io/ggeffects/index.html
The ggeffects package computes estimated marginal means (predicted values) for the response, at the margin of specific values or levels from certain model terms, i.e. it generates predictions by a model by holding the non-focal variables constant and varying the focal variable(s).
First try this for the models that don’t include random effects (fmod13 and fmod14)
#predict(fmod13,ForagePreen_wide,se.fit=TRUE) #couold use predict fxn but ggpredict makes it easier to graph:
#just one variable in model
fmod13_predict <- ggpredict(fmod13, terms = "feedingperhour_s")
ggplot(fmod13_predict, aes(x, predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)
plot(ggpredict(fmod13, terms = "feedingperhour_s"))
fmod13_predict2 <- ggpredict(fmod13, terms = c("feedingperhour_s","flocksize_nearest_s")) #both model terms
ggplot(fmod13_predict2, aes(x = x, y = predicted,color=group)) +
geom_line() 
plot(fmod13_predict2)
#with interaction:
#note the interaction is not significant in this model
#just one variable in model
fmod14_predict <- ggpredict(fmod14, terms = "feedingperhour_s")
ggplot(fmod14_predict, aes(x, predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)
plot(ggpredict(fmod14, terms = "feedingperhour_s"))
fmod14_predict2 <- ggpredict(fmod14, terms = c("feedingperhour_s","flocksize_nearest_s")) #both model terms
ggplot(fmod14_predict2, aes(x = x, y = predicted,color=group)) +
geom_line() 
plot(fmod14_predict2)
The models look very similar for both fmod13 and fmod14. When you look at the model summary for fmod14, the interaction is not significant. The interaction adds another parameter to the model (this interaction is a ‘pretending variable’). So they’re both supported but fmod14 has more parameters.
Now try this for models with site as random effect (fmod17 and fmod18)
#just one variable in model
fmod17_predict <- ggpredict(fmod17, terms = "feedingperhour_s")
ggplot(fmod17_predict, aes(x, predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)
plot(ggpredict(fmod17, terms = "feedingperhour_s"))
fmod17_predict2 <- ggpredict(fmod17, terms = c("feedingperhour_s","flocksize_nearest_s")) #both model terms
ggplot(fmod17_predict2, aes(x = x, y = predicted,color=group)) +
geom_line() 
plot(fmod17_predict2)
#with interaction:
#note the interaction is not significant in this model
#just one variable in model
fmod18_predict <- ggpredict(fmod18, terms = "feedingperhour_s")
ggplot(fmod14_predict, aes(x, predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)
plot(ggpredict(fmod18, terms = "feedingperhour_s")) 
plot(ggpredict(fmod18, terms = "flocksize_nearest_s"))
fmod18_predict2 <- ggpredict(fmod18, terms = c("feedingperhour_s","flocksize_nearest_s")) #both model terms
ggplot(fmod18_predict2, aes(x = x, y = predicted,color=group)) +
geom_line() 
plot(fmod18_predict2)
For the paper, showing model outputs for fmod17 (flock size + provisioning) might be useful. Choosing fmod17 because it doesn’t include the interaction in fmod18 which isn’t significant and does include site as random effect (as opposed to fmod13)
Make a 2 panel figure showing Foraging on Y axis and 1) provisioning and 2) flock size on X axis.
head(fmod17_predict) #data frame for provisioning## # Predicted counts of Foraging
##
## feedingperhour_s | Predicted | 95% CI
## -----------------------------------------------
## -1.26 | 252.83 | [192.44, 332.17]
## -0.75 | 229.93 | [185.30, 285.32]
## -0.41 | 216.15 | [178.51, 261.71]
## -0.33 | 212.78 | [176.44, 256.62]
## 0.49 | 182.75 | [148.85, 224.38]
## 0.54 | 181.09 | [146.90, 223.23]
##
## Adjusted for:
## * flocksize_nearest_s = 0.00
## * duration = 533.29
## * site = NA (population-level)fmod17_predict3 <- ggpredict(fmod17, terms = "flocksize_nearest_s")
head(fmod17_predict3)## # Predicted counts of Foraging
##
## flocksize_nearest_s | Predicted | 95% CI
## --------------------------------------------------
## -1.00 | 255.35 | [195.34, 333.80]
## -0.80 | 243.22 | [191.35, 309.15]
## -0.60 | 231.67 | [186.71, 287.46]
## -0.40 | 220.67 | [181.19, 268.74]
## -0.20 | 210.18 | [174.61, 253.01]
## 0.00 | 200.20 | [166.85, 240.21]
##
## Adjusted for:
## * feedingperhour_s = 0.00
## * duration = 533.29
## * site = NA (population-level)#merge these two dfs and then facet
#or just make each one pretty on it's own and paste together
p1<-ggplot(fmod17_predict, aes(x, predicted)) +
geom_line() +
geom_ribbon(data=fmod17_predict,aes(ymin = conf.low, ymax = conf.high), alpha = .1)+xlab("Provisioning Frequency (standardized)") + ylab("Predicted Foraging Time")+ggtitle("flock_size + provisioning + (1 | site)")+theme_classic()
p2<-ggplot(fmod17_predict3, aes(x, predicted)) +
geom_line() +
geom_ribbon(data=fmod17_predict3,aes(ymin = conf.low, ymax = conf.high), alpha = .1)+xlab("Flock Size (standardized)") + ylab(NULL)+theme_classic()
p12<-grid.arrange(p1,p2,ncol=2)
ggsave("outputs/foragingmodeloutput.png",plot=p12)## Saving 7 x 5 in image#make a plot that includes raw data and shows model fit:
#data = ForagePreen_wide
#this is Figure 3 in Manuscript
p1_points<-ggplot(fmod17_predict, aes(x, predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2,fill='#440154')+xlab("Provisioning frequency (standardized)") + ylab("Foraging time (seconds)")+ggtitle("flock_size + provisioning + (1 | site)")+theme_classic()+geom_point(data=ForagePreen_wide, aes(x=feedingperhour_s, y=Foraging),color='#ABABAB')
p2_points<-ggplot(fmod17_predict3, aes(x, predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2,fill='#440154')+xlab("Flock size (standardized)") + ylab("Foraging time (seconds)")+theme_classic()+geom_point(data=ForagePreen_wide, aes(x=flocksize_nearest_s, y=Foraging),color='#ABABAB')
fig3<-plot_grid(p1_points,p2_points, labels = c('A','B'))
#p12_points<-grid.arrange(p1_points,p2_points,ncol=2)
ggsave2("outputs/foragingmodeloutput_points.png",plot=fig3)## Saving 7 x 5 in imageMake the same plots but without standardizing the predictor variables.
fmod17_unscaled <-glmmTMB(Foraging~flocksize_nearest+feedingperhour+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(fmod17)## Family: nbinom2 ( log )
## Formula:
## Foraging ~ flocksize_nearest_s + feedingperhour_s + (1 | site) +
## offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1567.2 1585.3 -777.6 1555.2 144
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 2.283e-09 4.779e-05
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.97975 0.09296 -10.539 <2e-16 ***
## flocksize_nearest_s -0.24333 0.10301 -2.362 0.0182 *
## feedingperhour_s -0.18472 0.08670 -2.131 0.0331 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0330 0.1878 -5.499 3.82e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(fmod17_unscaled) #p values look the same## Family: nbinom2 ( log )
## Formula: Foraging ~ flocksize_nearest + feedingperhour + (1 | site) +
## offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1567.2 1585.3 -777.6 1555.2 144
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 2.52e-09 5.02e-05
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.534394 0.160263 -3.334 0.000855 ***
## flocksize_nearest -0.004545 0.001918 -2.369 0.017822 *
## feedingperhour -0.344850 0.161854 -2.131 0.033119 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0330 0.1878 -5.499 3.82e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fmod17_predict_unscaled <- ggpredict(fmod17_unscaled, terms = "feedingperhour")
fmod17_predict2_unscaled <- ggpredict(fmod17_unscaled, terms = "flocksize_nearest")
p1_points_unscaled<-ggplot(fmod17_predict_unscaled, aes(x, predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2,fill='#440154')+xlab("Provisioning Frequency") + ylab("Foraging Time (seconds)")+ggtitle("A).\n flock_size + provisioning + (1 | site)")+theme_classic()+geom_point(data=ForagePreen_wide, aes(x=feedingperhour, y=Foraging),color='#ABABAB')
p2_points_unscaled<-ggplot(fmod17_predict2_unscaled, aes(x, predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2,fill='#440154')+xlab("Flock Size") + ggtitle("B).")+ylab("Foraging Time (seconds)")+theme_classic()+geom_point(data=ForagePreen_wide, aes(x=flocksize_nearest, y=Foraging),color='#ABABAB')
p12_points_unscaled<-grid.arrange(p1_points_unscaled,p2_points_unscaled,ncol=2)
ggsave("foragingmodeloutput_points_unscaled.png",plot=p12_points_unscaled)## Saving 7 x 5 in imageDiagnostics using Dharma
Testing the fit of fmod17 and fmod18- the two best fitting models that include site as a random effect.
First test the fit of fmod17
## Family: nbinom2 ( log )
## Formula:
## Foraging ~ flocksize_nearest_s + feedingperhour_s + (1 | site) +
## offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1567.2 1585.3 -777.6 1555.2 144
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 2.283e-09 4.779e-05
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.97975 0.09296 -10.539 <2e-16 ***
## flocksize_nearest_s -0.24333 0.10301 -2.362 0.0182 *
## feedingperhour_s -0.18472 0.08670 -2.131 0.0331 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0330 0.1878 -5.499 3.82e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.65756, p-value = 0.148
## alternative hypothesis: two.sided
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0019, p-value = 1
## alternative hypothesis: two.sidedResiduals look fine. Plotting residuals against the two predictors suggests that there might be some issues for one of them.
Now test fit for fmod18
## Family: nbinom2 ( log )
## Formula:
## Foraging ~ flocksize_nearest_s * feedingperhour_s + (1 | site) +
## offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1569.1 1590.2 -777.6 1555.1 143
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 2.117e-09 4.601e-05
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.07
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.97710 0.09345 -10.455 <2e-16 ***
## flocksize_nearest_s -0.24949 0.10501 -2.376 0.0175 *
## feedingperhour_s -0.20032 0.09952 -2.013 0.0441 *
## flocksize_nearest_s:feedingperhour_s -0.05511 0.17401 -0.317 0.7515
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0333 0.1879 -5.499 3.81e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.66608, p-value = 0.182
## alternative hypothesis: two.sided
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0021, p-value = 1
## alternative hypothesis: two.sidedThe residuals look pretty okay for this model too but there are more problems when I plot against specific predictors.
Tables for Publication
Cand.models <- list("null + (1|site)" =fmod0.2,
"provisioning + (1|site)" = fmod15,
"flock size + (1|site)" = fmod16,
"provisioning + flock size + (1|site)" = fmod17,
"provisioning * flock size + (1|site)" = fmod18)
selectionTable <- aictab(cand.set = Cand.models, second.ord = FALSE) #uses AIC
row.names(selectionTable) <- NULL
kable(selectionTable,col.names = c("Model","K","AIC" , "$\\Delta$ AIC","ModelLik","AICWt", "LogLik","Cum.Wt"),caption = "<b>Foraging Candidate Models<b>") %>%
kable_styling(latex_options = "striped") %>%
kable_classic( lightable_options = "striped",full_width = F, html_font = "Cambria")| Model | K | AIC | \(\Delta\) AIC | ModelLik | AICWt | LogLik | Cum.Wt |
|---|---|---|---|---|---|---|---|
| provisioning + flock size + (1|site) | 6 | 1567.219 | 0.000000 | 1.0000000 | 0.4603299 | -777.6095 | 0.4603299 |
| provisioning * flock size + (1|site) | 7 | 1569.121 | 1.901545 | 0.3864425 | 0.1778910 | -777.5603 | 0.6382209 |
| flock size + (1|site) | 5 | 1569.594 | 2.374851 | 0.3050054 | 0.1404031 | -779.7969 | 0.7786240 |
| provisioning + (1|site) | 5 | 1569.842 | 2.623076 | 0.2694053 | 0.1240153 | -779.9211 | 0.9026393 |
| null + (1|site) | 4 | 1570.326 | 3.107042 | 0.2115019 | 0.0973607 | -781.1630 | 1.0000000 |
sf17<-summary(fmod17)
sumfmod17<-as.data.frame(sf17$coefficients[1])
kable(sumfmod17,col.names = c("Estimate","Std. error","z value","p-value"),caption = "<b>fmod17 summary<b>") %>%
kable_styling(latex_options = "striped") %>%
kable_classic( lightable_options = "striped",full_width = F, html_font = "Cambria")| Estimate | Std. error | z value | p-value | |
|---|---|---|---|---|
| (Intercept) | -0.9797530 | 0.0929608 | -10.539422 | 0.0000000 |
| flocksize_nearest_s | -0.2433306 | 0.1030108 | -2.362186 | 0.0181675 |
| feedingperhour_s | -0.1847193 | 0.0866983 | -2.130599 | 0.0331222 |
sf18<-summary(fmod18)
sumfmod18<-as.data.frame(sf18$coefficients[1])
kable(sumfmod18,col.names = c("Estimate","Std. error","z value","p-value"),caption = "<b>fmod18 summary<b>") %>%
kable_styling(latex_options = "striped") %>%
kable_classic( lightable_options = "striped",full_width = F, html_font = "Cambria")| Estimate | Std. error | z value | p-value | |
|---|---|---|---|---|
| (Intercept) | -0.9770989 | 0.0934536 | -10.4554490 | 0.0000000 |
| flocksize_nearest_s | -0.2494926 | 0.1050081 | -2.3759363 | 0.0175045 |
| feedingperhour_s | -0.2003245 | 0.0995204 | -2.0128988 | 0.0441253 |
| flocksize_nearest_s:feedingperhour_s | -0.0551059 | 0.1740110 | -0.3166803 | 0.7514862 |
#try to merge these two summary tables
sumfmod17$Model <- "fmod17"
sumfmod18$Model <- "fmod18"
sumfmod17$Predictors <- rownames(sumfmod17) #convert rownames to column
sumfmod18$Predictors <- rownames(sumfmod18)#convert rownames to column
fmods <- rbind(sumfmod17, sumfmod18)
rownames(fmods)<-NULL #remove rownames
fmods <- fmods %>%
dplyr::select(Predictors, everything())
fmods$Model<-NULL
fmods <- fmods %>%
rename("pval" = "cond.Pr...z..") #renaming column
#fmods$pval<-
# dplyr::select(fmods,pval)%>% mutate_if(is.numeric, ~ case_when(. <0.00000001 ~ round(., 30), TRUE ~ round(.,5)))
#fmods <- fmods %>%
# rename("pval" = "cond.Pr...z..") #renaming column
#fmods <- fmods %>%
# mutate(pval = case_when(
# pval < 0.05 ~ paste0(as.character(pval), "*"),
# TRUE ~ paste0(pval)
# ))
kable(fmods,col.names = c("Response/Predictors","Estimate","Std. error","z-value","p-value"),digits = c(5,5,5,5,20)) %>%
kable_styling(latex_options = "striped") %>%
kable_classic( lightable_options = "striped",full_width = F, html_font = "Cambria")%>%
pack_rows("Foraging", 1, 3) %>%
pack_rows("Foraging", 4, 7)| Response/Predictors | Estimate | Std. error | z-value | p-value |
|---|---|---|---|---|
| Foraging | ||||
| (Intercept) | -0.97975 | 0.09296 | -10.53942 | 0.00000000 |
| flocksize_nearest_s | -0.24333 | 0.10301 | -2.36219 | 0.01816752 |
| feedingperhour_s | -0.18472 | 0.08670 | -2.13060 | 0.03312218 |
| Foraging | ||||
| (Intercept) | -0.97710 | 0.09345 | -10.45545 | 0.00000000 |
| flocksize_nearest_s | -0.24949 | 0.10501 | -2.37594 | 0.01750448 |
| feedingperhour_s | -0.20032 | 0.09952 | -2.01290 | 0.04412528 |
| flocksize_nearest_s:feedingperhour_s | -0.05511 | 0.17401 | -0.31668 | 0.75148618 |
Grooming models:
Main Q:
Does flock size, provisioning, both, or neither influence foraging time?
Grooming includes preening or bathing.
Exploratory Visualizations
Model selection:
Count data: Does poisson or negative binomial fit better? Build global model (two predictors + interactions) and compare AIC.
pm1_preen = glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s + offset(log(duration)), ForagePreen_wide, family=poisson)
nbm1_preen = glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), ForagePreen_wide, family=nbinom2)
zipm1_preen = glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~1, ForagePreen_wide, family=poisson)
zinbm1_preen = glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~1, ForagePreen_wide, family=nbinom2)
AICtab(pm1_preen,nbm1_preen,zipm1_preen,zinbm1_preen) #zero inflated is best here, continue with this model## dAIC df
## zinbm1_preen 0.0 6
## nbm1_preen 31.7 5
## zipm1_preen 14544.6 5
## pm1_preen 24584.0 4Zero inflated negative binomial model is the best one so continue with that model structure.
##No random effects:
#intercept only
pmod0.3 <-glmmTMB(Preening~1+ offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(pmod0.3)## Family: nbinom2 ( log )
## Formula: Preening ~ 1 + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1545.9 1554.9 -769.9 1539.9 147
##
##
## Dispersion parameter for nbinom2 family (): 0.763
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.158 0.108 -10.72 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1912 0.2037 -5.849 4.95e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning
pmod11 <-glmmTMB(Preening~feedingperhour_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(pmod11)## Family: nbinom2 ( log )
## Formula: Preening ~ feedingperhour_s + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1547.2 1559.3 -769.6 1539.2 146
##
##
## Dispersion parameter for nbinom2 family (): 0.766
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.16616 0.10803 -10.795 <2e-16 ***
## feedingperhour_s 0.09753 0.12266 0.795 0.427
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1919 0.2038 -5.848 4.98e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#flock size
pmod12 <-glmmTMB(Preening~flocksize_nearest_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(pmod12)## Family: nbinom2 ( log )
## Formula: Preening ~ flocksize_nearest_s + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1545.0 1557.0 -768.5 1537.0 146
##
##
## Dispersion parameter for nbinom2 family (): 0.783
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1765 0.1067 -11.022 <2e-16 ***
## flocksize_nearest_s 0.1638 0.1016 1.613 0.107
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1870 0.2027 -5.856 4.74e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning + flock size
pmod13 <-glmmTMB(Preening~flocksize_nearest_s+feedingperhour_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(pmod13)## Family: nbinom2 ( log )
## Formula:
## Preening ~ flocksize_nearest_s + feedingperhour_s + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1546.6 1561.6 -768.3 1536.6 145
##
##
## Dispersion parameter for nbinom2 family (): 0.784
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.18239 0.10677 -11.074 <2e-16 ***
## flocksize_nearest_s 0.15642 0.10063 1.554 0.120
## feedingperhour_s 0.07797 0.11909 0.655 0.513
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1875 0.2028 -5.855 4.76e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning x flock size
pmod14 <-glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(pmod14)## Family: nbinom2 ( log )
## Formula:
## Preening ~ flocksize_nearest_s * feedingperhour_s + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1548.6 1566.6 -768.3 1536.6 144
##
##
## Dispersion parameter for nbinom2 family (): 0.784
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.18402 0.10912 -10.851 <2e-16 ***
## flocksize_nearest_s 0.15785 0.10248 1.540 0.123
## feedingperhour_s 0.08129 0.12786 0.636 0.525
## flocksize_nearest_s:feedingperhour_s 0.01284 0.17902 0.072 0.943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1877 0.2029 -5.855 4.77e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AICtab(pmod0.3,pmod11,pmod12,pmod13,pmod14)## dAIC df
## pmod12 0.0 4
## pmod0.3 0.9 3
## pmod13 1.6 5
## pmod11 2.2 4
## pmod14 3.6 6##Same models but with site as random effect
#intercept only
pmod0.2 <-glmmTMB(Preening~1+ (1|site) + offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(pmod0.2)## Family: nbinom2 ( log )
## Formula: Preening ~ 1 + (1 | site) + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1545.8 1557.8 -768.9 1537.8 146
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.07115 0.2667
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 0.802
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2169 0.1624 -7.491 6.84e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1857 0.2025 -5.856 4.73e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning
pmod15 <-glmmTMB(Preening~feedingperhour_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(pmod15)## Family: nbinom2 ( log )
## Formula:
## Preening ~ feedingperhour_s + (1 | site) + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1547.7 1562.8 -768.9 1537.7 145
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.08273 0.2876
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 0.805
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.22063 0.17031 -7.167 7.65e-13 ***
## feedingperhour_s -0.03225 0.17501 -0.184 0.854
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1851 0.2023 -5.857 4.72e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#flock size
pmod16 <-glmmTMB(Preening~flocksize_nearest_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(pmod16)## Family: nbinom2 ( log )
## Formula:
## Preening ~ flocksize_nearest_s + (1 | site) + offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1546.6 1561.7 -768.3 1536.6 145
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.03272 0.1809
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 0.799
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2026 0.1384 -8.686 <2e-16 ***
## flocksize_nearest_s 0.1329 0.1219 1.091 0.275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1848 0.2022 -5.858 4.68e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning + flock size
pmod17 <-glmmTMB(Preening~flocksize_nearest_s+feedingperhour_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(pmod17)## Family: nbinom2 ( log )
## Formula:
## Preening ~ flocksize_nearest_s + feedingperhour_s + (1 | site) +
## offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1548.5 1566.6 -768.2 1536.5 144
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01502 0.1226
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 0.793
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.19393 0.12748 -9.366 <2e-16 ***
## flocksize_nearest_s 0.14585 0.11664 1.250 0.211
## feedingperhour_s 0.05801 0.15247 0.380 0.704
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1861 0.2026 -5.855 4.77e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning x flock size
pmod18 <-glmmTMB(Preening~flocksize_nearest_s*feedingperhour_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=ForagePreen_wide)
summary(pmod18)## Family: nbinom2 ( log )
## Formula:
## Preening ~ flocksize_nearest_s * feedingperhour_s + (1 | site) +
## offset(log(duration))
## Zero inflation: ~1
## Data: ForagePreen_wide
##
## AIC BIC logLik deviance df.resid
## 1550.5 1571.6 -768.2 1536.5 143
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01476 0.1215
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 0.792
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.194963 0.128629 -9.290 <2e-16 ***
## flocksize_nearest_s 0.147157 0.118857 1.238 0.216
## feedingperhour_s 0.060848 0.161153 0.378 0.706
## flocksize_nearest_s:feedingperhour_s 0.009679 0.179544 0.054 0.957
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1862 0.2026 -5.855 4.78e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AICtab(pmod0.2,pmod15,pmod16,pmod17,pmod18)## dAIC df
## pmod0.2 0.0 4
## pmod16 0.8 5
## pmod15 2.0 5
## pmod17 2.7 6
## pmod18 4.7 7The top performing models are within 2 AIC values of the intercept only null model.
Model diagnostics
Not needed because the intercept only model is the best fitting model?
Tables for Publication
Cand.models_groom <- list("null + (1|site)" =pmod0.2,
"provisioning + (1|site)" = pmod15,
"flock size + (1|site)" = pmod16,
"provisioning + flock size + (1|site)" = pmod17,
"provisioning * flock size + (1|site)" = pmod18)
selectionTable_groom <- aictab(cand.set = Cand.models_groom, second.ord = FALSE) #uses AIC
row.names(selectionTable_groom) <- NULL
kable(selectionTable_groom,col.names = c("Model","K","AIC" , "$\\Delta$ AIC","ModelLik","AICWt", "LogLik","Cum.Wt"),caption = "<b>Grooming Candidate Models<b>") %>%
kable_styling(latex_options = "striped") %>%
kable_classic( lightable_options = "striped",full_width = F, html_font = "Cambria")| Model | K | AIC | \(\Delta\) AIC | ModelLik | AICWt | LogLik | Cum.Wt |
|---|---|---|---|---|---|---|---|
| null + (1|site) | 4 | 1545.776 | 0.0000000 | 1.0000000 | 0.4200894 | -768.8879 | 0.4200894 |
| flock size + (1|site) | 5 | 1546.623 | 0.8474312 | 0.6546100 | 0.2749948 | -768.3116 | 0.6950842 |
| provisioning + (1|site) | 5 | 1547.742 | 1.9657806 | 0.3742279 | 0.1572092 | -768.8708 | 0.8522934 |
| provisioning + flock size + (1|site) | 6 | 1548.494 | 2.7177867 | 0.2569450 | 0.1079399 | -768.2468 | 0.9602332 |
| provisioning * flock size + (1|site) | 7 | 1550.491 | 4.7148714 | 0.0946627 | 0.0397668 | -768.2454 | 1.0000000 |
Vigilance models:
Main Q:
Does flock size, provisioning, both, or neither influence vigilance time?
Exploratory Visualizations
Model selection:
Count data: Does poisson or negative binomial fit better? Build global model (two predictors + interactions) and compare AIC.
pm1_vig = glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s + offset(log(duration)), Vig_wide, family=poisson)
nbm1_vig = glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), Vig_wide, family=nbinom2)
zipm1_vig = glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~1, Vig_wide, family=poisson)
zinbm1_vig = glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s +offset(log(duration)), zi=~1, Vig_wide, family=nbinom2)
AICtab(pm1_vig,nbm1_vig,zipm1_vig,zinbm1_vig) #zero inflated is best here, continue with this model## dAIC df
## zinbm1_vig 0.0 6
## nbm1_vig 14.0 5
## zipm1_vig 17292.7 5
## pm1_vig 19778.4 4Zero inflated negative binomial model is the best one so continue with that model structure.
##No random effects:
#intercept only
vmod0.3 <-glmmTMB(Vigilant~1+ offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=Vig_wide)
summary(vmod0.3)## Family: nbinom2 ( log )
## Formula: Vigilant ~ 1 + offset(log(duration))
## Zero inflation: ~1
## Data: Vig_wide
##
## AIC BIC logLik deviance df.resid
## 1795.1 1804.1 -894.5 1789.1 147
##
##
## Dispersion parameter for nbinom2 family (): 1.04
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1186 0.0831 -13.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7527 0.3674 -7.493 6.72e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning
vmod11 <-glmmTMB(Vigilant~feedingperhour_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=Vig_wide)
summary(vmod11)## Family: nbinom2 ( log )
## Formula: Vigilant ~ feedingperhour_s + offset(log(duration))
## Zero inflation: ~1
## Data: Vig_wide
##
## AIC BIC logLik deviance df.resid
## 1795.2 1807.3 -893.6 1787.2 146
##
##
## Dispersion parameter for nbinom2 family (): 1.05
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.12883 0.08280 -13.634 <2e-16 ***
## feedingperhour_s 0.11559 0.08451 1.368 0.171
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7571 0.3691 -7.471 7.98e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#flock size
vmod12 <-glmmTMB(Vigilant~flocksize_nearest_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=Vig_wide)
summary(vmod12)## Family: nbinom2 ( log )
## Formula: Vigilant ~ flocksize_nearest_s + offset(log(duration))
## Zero inflation: ~1
## Data: Vig_wide
##
## AIC BIC logLik deviance df.resid
## 1796.7 1808.7 -894.3 1788.7 146
##
##
## Dispersion parameter for nbinom2 family (): 1.05
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.11978 0.08298 -13.494 <2e-16 ***
## flocksize_nearest_s -0.05246 0.08319 -0.631 0.528
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7504 0.3665 -7.504 6.21e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning + flock size
vmod13 <-glmmTMB(Vigilant~flocksize_nearest_s+feedingperhour_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=Vig_wide)
summary(vmod13)## Family: nbinom2 ( log )
## Formula:
## Vigilant ~ flocksize_nearest_s + feedingperhour_s + offset(log(duration))
## Zero inflation: ~1
## Data: Vig_wide
##
## AIC BIC logLik deviance df.resid
## 1796.6 1811.7 -893.3 1786.6 145
##
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.13097 0.08263 -13.688 <2e-16 ***
## flocksize_nearest_s -0.06498 0.08119 -0.800 0.423
## feedingperhour_s 0.12392 0.08560 1.448 0.148
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7546 0.3682 -7.482 7.33e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning x flock size
vmod14 <-glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=Vig_wide)
summary(vmod14)## Family: nbinom2 ( log )
## Formula:
## Vigilant ~ flocksize_nearest_s * feedingperhour_s + offset(log(duration))
## Zero inflation: ~1
## Data: Vig_wide
##
## AIC BIC logLik deviance df.resid
## 1798.5 1816.6 -893.2 1786.5 144
##
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.13733 0.08423 -13.503 <2e-16 ***
## flocksize_nearest_s -0.05786 0.08320 -0.695 0.487
## feedingperhour_s 0.13835 0.09467 1.461 0.144
## flocksize_nearest_s:feedingperhour_s 0.05327 0.14684 0.363 0.717
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7558 0.3687 -7.474 7.75e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AICtab(vmod0.3,vmod11,vmod12,vmod13,vmod14)## dAIC df
## vmod0.3 0.0 3
## vmod11 0.2 4
## vmod13 1.5 5
## vmod12 1.6 4
## vmod14 3.4 6##Same models but with site as random effect
#intercept only
vmod0.2 <-glmmTMB(Vigilant~1+ (1|site) + offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=Vig_wide)
summary(vmod0.2)## Family: nbinom2 ( log )
## Formula: Vigilant ~ 1 + (1 | site) + offset(log(duration))
## Zero inflation: ~1
## Data: Vig_wide
##
## AIC BIC logLik deviance df.resid
## 1796.7 1808.8 -894.4 1788.7 146
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01613 0.127
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.1253 0.1008 -11.16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7493 0.3662 -7.508 5.99e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning
vmod15 <-glmmTMB(Vigilant~feedingperhour_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=Vig_wide)
summary(vmod15)## Family: nbinom2 ( log )
## Formula:
## Vigilant ~ feedingperhour_s + (1 | site) + offset(log(duration))
## Zero inflation: ~1
## Data: Vig_wide
##
## AIC BIC logLik deviance df.resid
## 1796.9 1812.0 -893.5 1786.9 145
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.0167 0.1292
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.07
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.13354 0.10079 -11.247 <2e-16 ***
## feedingperhour_s 0.12838 0.09711 1.322 0.186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7515 0.3671 -7.495 6.61e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#flock size
vmod16 <-glmmTMB(Vigilant~flocksize_nearest_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=Vig_wide)
summary(vmod16)## Family: nbinom2 ( log )
## Formula:
## Vigilant ~ flocksize_nearest_s + (1 | site) + offset(log(duration))
## Zero inflation: ~1
## Data: Vig_wide
##
## AIC BIC logLik deviance df.resid
## 1798.4 1813.4 -894.2 1788.4 145
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01607 0.1268
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.12738 0.10085 -11.18 <2e-16 ***
## flocksize_nearest_s -0.05547 0.09087 -0.61 0.542
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7480 0.3657 -7.514 5.73e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning + flock size
vmod17 <-glmmTMB(Vigilant~flocksize_nearest_s+feedingperhour_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=Vig_wide)
summary(vmod17)## Family: nbinom2 ( log )
## Formula:
## Vigilant ~ flocksize_nearest_s + feedingperhour_s + (1 | site) +
## offset(log(duration))
## Zero inflation: ~1
## Data: Vig_wide
##
## AIC BIC logLik deviance df.resid
## 1798.5 1816.6 -893.2 1786.5 144
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.01066 0.1032
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.07
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.13444 0.09477 -11.971 <2e-16 ***
## flocksize_nearest_s -0.05885 0.08814 -0.668 0.504
## feedingperhour_s 0.12762 0.09276 1.376 0.169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7515 0.3671 -7.495 6.64e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#provisioning x flock size
vmod18 <-glmmTMB(Vigilant~flocksize_nearest_s*feedingperhour_s+(1|site)+offset(log(duration)),
family=nbinom2,
ziformula=~1,
data=Vig_wide)
summary(vmod18)## Family: nbinom2 ( log )
## Formula:
## Vigilant ~ flocksize_nearest_s * feedingperhour_s + (1 | site) +
## offset(log(duration))
## Zero inflation: ~1
## Data: Vig_wide
##
## AIC BIC logLik deviance df.resid
## 1800.4 1821.5 -893.2 1786.4 143
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.007085 0.08417
## Number of obs: 150, groups: site, 5
##
## Dispersion parameter for nbinom2 family (): 1.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.13790 0.09207 -12.359 <2e-16 ***
## flocksize_nearest_s -0.05544 0.08769 -0.632 0.527
## feedingperhour_s 0.13647 0.09856 1.385 0.166
## flocksize_nearest_s:feedingperhour_s 0.03804 0.16156 0.235 0.814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7533 0.3679 -7.483 7.24e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AICtab(vmod0.2,vmod15,vmod16,vmod17,vmod18)## dAIC df
## vmod0.2 0.0 4
## vmod15 0.2 5
## vmod16 1.6 5
## vmod17 1.8 6
## vmod18 3.7 7The top performing models are within 2 AIC values of the intercept only null model.
Model diagnostics
Not needed because the intercept only model is the best fitting model?
Tables for Publication
Cand.models_vig <- list("null + (1|site)" =vmod0.2,
"provisioning + (1|site)" = vmod15,
"flock size + (1|site)" = vmod16,
"provisioning + flock size + (1|site)" = vmod17,
"provisioning * flock size + (1|site)" = vmod18)
selectionTable_vig <- aictab(cand.set = Cand.models_vig, second.ord = FALSE) #uses AIC
row.names(selectionTable_vig) <- NULL
kable(selectionTable_vig,col.names = c("Model","K","AIC" , "$\\Delta$ AIC","ModelLik","AICWt", "LogLik","Cum.Wt"),caption = "<b>Vigilance Candidate Models<b>") %>%
kable_styling(latex_options = "striped") %>%
kable_classic( lightable_options = "striped",full_width = F, html_font = "Cambria")| Model | K | AIC | \(\Delta\) AIC | ModelLik | AICWt | LogLik | Cum.Wt |
|---|---|---|---|---|---|---|---|
| null + (1|site) | 4 | 1796.749 | 0.0000000 | 1.0000000 | 0.3405399 | -894.3743 | 0.3405399 |
| provisioning + (1|site) | 5 | 1796.914 | 0.1650839 | 0.9207728 | 0.3135599 | -893.4569 | 0.6540997 |
| flock size + (1|site) | 5 | 1798.384 | 1.6356869 | 0.4413825 | 0.1503083 | -894.1922 | 0.8044081 |
| provisioning + flock size + (1|site) | 6 | 1798.500 | 1.7508891 | 0.4166767 | 0.1418950 | -893.2498 | 0.9463031 |
| provisioning * flock size + (1|site) | 7 | 1800.443 | 3.6943541 | 0.1576817 | 0.0536969 | -893.2215 | 1.0000000 |